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Abstract. In this work we present a novel framework for the computation of 
finite dimensional invariant sets of infinite dimensional dynamical systems. It 
extends a classical subdivision technique [5] for the computation of such objects 
of finite dimensional systems to the infinite dimensional case by utilizing results 
on embedding techniques for infinite dimensional systems. We show how to 
implement this approach for the analysis of delay differential equations and 
illustrate the feasibility of our implementation by computing invariant sets for 
three different delay differential equations. 
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1. Introduction 

Over the last two decades so-called set oriented numerical methods have been 
developed in the context of the numerical treatment of dynamical systems (e. g. 
[5, 6, 11, 13]). The basic idea is to cover the objects of interest - for instance 
invariant sets or invariant measures - by outer approximations which are created 
via multilevel subdivision techniques. These techniques have been used in several 
different application areas such as molecular dynamics ([23]), astrodynamics ([7]) 
or ocean dynamics ([12]). 

So far the applicability of the subdivision scheme is restricted to finite dimen¬ 
sional dynamical systems, i. e. ordinary differential equations or finite dimensional 
discrete dynamical systems. In this paper we extend this technique to the infinite 
dimensional context. More precisely, we develop a set oriented numerical technique 
which allows us to compute low-dimensional invariant sets for infinite dimensional 
dynamical systems. Rather than using a straightforward approach based on an 
appropriate combination of Galerkin expansions and subdivision steps we follow a 
completely novel path and utilize embedding results in our numerical treatment. 

The first result on embeddings in the dynamical systems context is the celebrated 
Tokens Embedding Theorem [25]. Takens has shown that an invariant set - in his 
context this has to be a compact manifold of dimension d - can generically be 
reconstructed using the so-called observation map which consists of observations of 
the dynamical behavior at an appropriate number (at least 2d + 1) of consecutive 
snapshots in time. This result has been extended by Sauer et al. in [22] to the context 
of compact invariant sets of box counting dimension d. There it has been shown 
that the same observation map can be used for the reconstruction of the invariant 
set as long as more than 2d consecutive snapshots in time are used. Moreover in 
this work the notion of “genericity” has been replaced by the more intuitive notion 
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of “prevalence”. Finally, in [20] Robinson extended the results obtained in [22] to 
dynamical systems on infinite dimensional Banach spaces. It turns out that also 
here the same observation map can be used to reconstruct invariant sets of finite 
dimensional box counting dimension d. However, in addition to the dimension of 
the set another quantity comes into play, namely the thickness exponent a. Roughly 
speaking this exponent measures how well the invariant set can be approximated 
by finite dimensional subspaces of the underlying Banach space. The lower bound 
2d on the number of snapshots has to be replaced by 2(1 + a)d accordingly. 

We remark that there are several further extensions of Takens’ theorem. For 
instance in [24] forced systems are considered, and in [19] one can find a stochastic 
version of this result. 

Our results in this paper are based on Robinson’s embedding theorem. In fact, 
we will combine the reconstruction based on the observation map with the classical 
subdivision techniques developed in [5] . Assuming that a bound on the box-counting 
dimension and the thickness exponent of the invariant set are known we use the 
observation map and its inverse to define a dynamical system cp in the embedding 
space of dimension k > 2(1 -\-a)d. Then the subdivision scheme is applied to compute 
the reconstructed invariant set for cp. Observe that this way we can always perform 
the numerical approximation within a finite dimensional space of fixed dimension 
&, and this is in contrast to Galerkin based approaches where one would have to 
extend the expansions in order to improve the quality of the approximation. 

The numerical approach we are proposing is in principle applicable to arbitrary 
infinite dimensional dynamical systems. However, here we will restrict our attention 
to delay differential equations with constant delay in the numerical realization. The 
applications to e. g. partial differential equations will be done in future work. 

A detailed outline of the paper is as follows. In Section 2 we briefly summarize 
the infinite dimensional embedding theory introduced in [20]. In Section 3 we 
employ the main result of [20] for the construction of a numerical approach for 
the computation of compact finite dimensional attractors of infinite dimensional 
dynamical systems. First we construct a continuous dynamical system cp on the 
embedding space using a generalization of the well-known Tietze extension theorem 
[9, 1.5.3] which is due to Dugundji [8, Theorem 4.1]. Then we extend the results from 
[5] to the situation where the underlying dynamical system (p is just continuous (and 
not homoemorphic). A numerical realization for the computation of attractors of 
delay differential equations is given in Section 4. Finally, in Section 5, we illustrate 
the efficiency of our novel approach for three different delay differential equations. 

2. Infinite Dimensional Takens Embedding 

We start with a short review of the contents of [20]. We consider a dynamical 
system of the form 

(1) u j+1 = j = 0,1,..., 

where : Y -A Y is Lipschitz continuous on a Banach space Y. Moreover we 
assume that <f> has an invariant compact set A, that is 

<f>(A) = A. 

Later on we will additionally assume that A is a global attractor in the sense 
that it attracts all bounded sets within Y as t -A oo. 
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For the statement of the main result of [20] we need three particular notions: 
prevalence [22], upper box counting dimension and thickness exponent [15]. 

Definition 2.1 (Prevalence). A Borel subset S of a normed linear space V is 
prevalent if there is a finite dimensional subspace E of V (the ‘probe space’) such 
that for each v G V, v + e belongs to S for (Lebesgue) almost every e G E. 

Definition 2.2 (Upper box counting dimension). Let Y be a Banach space, and 
let X C Y be compact. For 5 > 0, denote by Ny{X , e) the minimal number of balls 
of radius 6 (in the norm of Y) necessary to cover the set X. Then 

(2) d(X]Y) — limsup ^°&Ny(X,£) p msU p —log Ny(X,e) 

£—^0 - log £ £ —^0 

denotes the upper box-counting dimension of X. 

Definition 2.3 (Thickness exponent). Let Y be a Banach space, and let X C Y 
be compact. For 6 > 0, denote by d(X,e) the minimal dimension of all finite 
dimensional subspaces V C Y such that every point of X lies within distance 5 of 
V; if no such V exists, d(X,e) = oc. Then 

cr(X,Y) := limsup —log £ d(X,e) 

£ —H) 

is called the thickness exponent of X in Y. 

From this definition one sees that essentially, cr(X,Y) captures how well X can 
be approximated from within finite dimensional subspaces of Y. In [17], as a more 
practical expression Kukavica and Robinson prove that 

cr(X, Y) = limsup—-— ^ ™ —-, i.e. ey(X, n) ~ n ~ 1 / cr ( x X) 
n—>-oo -log £y(X,n) 

where £y (A", n) is the minimum distance between X and any n-dimensional linear 
subspace of Y. 

These notions are essential in answering the question when a delay embedding 
technique applied to an invariant subset A <zY will work generically. More precisely, 
the results are as follows. 


Theorem 2.4 ([15]). Let Y be a Banach space and X C Y compact, with upper 
box counting dimension d{X\Y) =: d and thickness exponent cr(X,Y) =: a. Let 
N > 2d be an integer, and let a E M with 


0 < a < 


N -2d 
AU(l + a)‘ 


Then for a prevalent set of linear maps C :Y —>> R N there is C > 0 such that 


C • \\£{x-y)\\ a > \\x-y\\ for all x,y G X. 


Note that this implies that a prevalent set of linear maps Y —>• will be one-to- 
one on X. Using this theorem, the following result concerning the delay embedding 
technique can be proven. 

Theorem 2.5 ([20]). Suppose that the upper box counting dimension of A is d(A) = 
d, and that A has a thickness exponent a. Choose an integer k > 2(1 + cr)d and 
suppose further that the set A p of p-periodic points o/4> satisfies d(A p ) < p /(2 + 2 a) 





4 


M. DELLNITZ, M. HESSEL-VON MOLO A. ZIESSLER 


for p = 1 ,...,&. Then for a prevalent set of Lipschitz maps f : Y -A M 
observation map Dk[f , $] : T -A defined by 

(3) £>*[/, *](«) = (/(«), /(*(«)), • • •, /($ fe - 1 («))) T 

zs one-to-one on the invariant set A. 

Remark 1. Following an observation already made in [22, Remark 2.9], we note 
that this result may be generalized to the case where several independent observables 
are evaluated. More precisely, also for a prevalent set of Lipschitz maps /i,..., f q : 
Y -A M the observation map D k [fi ,..., f q \ : Y -a M fe , 

(4) « (/!(«), . . . , fq(u), f q (^~ 1 (u))) T 

is one-to-one on A, provided that k = J2i=i > 2(1 + a) • d and d(A p ) < p/{2 + 2cr) 
for p < max(fci,..., fc g ). 

3. Computation of Embedded Attractors via Subdivision 

3.1. Finite-dimensional Embeddings of Attractors. In this section we employ 
Theorem 2.5 in order to construct a method for the computation of compact, finite 
dimensional attractors of infinite dimensional dynamical systems on a Banach space 
Y. 

Let us denote by A^ the image of A cY under D [/, <f>], that is 

A k = D k [f,$](A), 

where D k is the map defined in Theorem 2.5 and / is chosen such that D k is 
one-to-one on A. 

We now develop a set oriented numerical technique for the approximation of the 
set A k . First, we define a dynamical system on for which A k is an invariant set, 
on which the dynamics is conjugate to that of <F on A. For this we define the map 
p : M. k -A M. k by 

(5) p = R o <F o E, 

where E : -A Y and R : Y -A are an embedding and a restriction, respectively, 

that satisfy 

(6) (E o R)(u) = u \/u G A and (RoE)(x)=x \/x E A k . 

More concretely, we define the map R = by the right-hand side of (3) (see 
Figure 1). 


E 



R=D k \f, O] 


Figure 1. Definition of the map p. 
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Remark 2. Although Dk is derived from what would commonly be termed an 
embedding theorem, we call R a restriction, as it maps from an infinite dimensional 
domain into a finite-dimensional set, and E an embedding, as it maps (or embeds) 
a finite-dimensional into an infinite dimensional space. 

The map E is obtained in two steps: Firstly, as R is one-to-one on A , the 
requirement that 

(7) R o E(x) = x for all x E A k 

uniquely defines a map E : Ak —> A. In a second step, we extend this map to 
a continuous map E : R k -A Y. To do this, we employ a generalization of the 
well-known Tietze extension theorem [9, 1.5.3] found by Dugundji [8, Theorem 4.1], 
stated here with notation adapted to our needs: 

Theorem 3.1. Let X be an arbitrary metric space and A C X be closed, V a locally 
convex linear space and p : A -A V be continuous. Then there is a continuous map 
P : X -A V with P\a = P such that P(X) is contained in the convex hull of p(A). 

Using this theorem in the situation introduced above, we obtain the following 

Proposition 1. There is a continuous map ip : R k -A R k satisfying 

(8) (p(R(u )) = R($(u)) for all u E A. 

Proof. By construction, the map R : Y -A R k given by (3) is continuous (even 
Lipschitz) and one-to-one. Thus, restricting R to A we obtain a bijective map 
R : A -a Ak. As A is assumed to be compact and Ak C R k is Hausdorff, R is a 
homeomorphism by a well-known theorem from elementary topology (see e. g. [26, 
Theorem 17.14]). Thus we obtain a continuous map E : A^ -A A as E = i? -1 . 

As Y is a normed linear space, it is locally convex. Thus we can apply Dugundji’s 
Theorem with X = R k , A = Ak, p = E and V = Y to see that there is a continuous 
map E : R k ^ Y with E\^ k = E. Finally, defining p through (5) we obtain that 
p is continuous as a composition of continous maps. □ 

Now we are in a position to approximate the embedded invariant set Ak via the 
corresponding dynamical system 

(9) x j+1 =<p(xj) j= 0,1,.... 

To this end, we employ a subdivision scheme as defined in [5]. 

3.2. Subdivision Scheme. We briefly review the classical subdivision procedure. 
Let Q C be a compact set. We define the global attractor relative to Q by 

( 10 ) A Q =f] P{Q). 

j> 0 

The subdivision procedure allows us to approximate this set. Roughly speaking, the 
idea of the algorithm is as follows. We start with a finite family of (large) compact 
subsets of R k which cover the domain in which we want to analyze the dynamical 
behavior. Then we subdivide each of these sets into smaller ones and throw away 
subsets which do not contain part of the relative global attractor. Continuing the 
process with the new collection of (smaller) sets it becomes intuitively clear that this 
should lead to a successively finer approximation of the relative global attractor. 
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Let us be more precise. The algorithm generates a sequence So, Si,... of finite 
collections of compact subsets of R k such that the diameter 

diam(S^) = max diam(S) 

Bet3i 

converges to zero for i -A oo. Given an initial collection So, we inductively obtain 
Bg from Bi -1 for i = 1,2,... in two steps. 

(1) Subdivision: Construct a new collection B>i such that 

oo U B = U B 

BeBe BeBi-! 

and 

(12) diam(S^) = diam(S^_i), 

where 0 <C Gnin ^ — ^max ^ L 

(2) Selection: Define the new collection Bi by 

(13) Bi = |s G B>i : 3S G S^ such that cp~ 1 (B) D B 0| . 

The first step guarantees that the collections Bn consist of successively finer sets 
for increasing £. In fact, by construction 

diam(S^) < 6^ ax diam(So) —^ 0 for £ -A oo. 

In the second step we remove each subset whose preimage does neither intersect 
itself nor any other subset in B>£. As we shall see, this step is responsible for the 
fact that the unions \J Be g e B approach the relative global attractor. 

Denote by Qa the collection of compact subsets obtained after t, subdivision steps, 
that is, 

Qe= |J B - 

BeBi 

Moreover let Bq be a finite collection of closed subsets with Q o = (J beb 0 B = Q. 
Then the main convergence result of [5] states that 

lim h (Aq, Qk) = 0, 

k—yoo 

where h(B,C) is the usual Hausdorff distance between two compact subsets B,C C 
M n . However, in that work the authors assume that cp is a homeomorphism and 
not just continuous, as in the situation here. For this reason, in the following we 
present a proof of convergence for continuous ip. 

3.3. Proof of Convergence. Essentially, we will be able to follow the structure 
of the proof in [5]. However, there are some technical differences, and we will need 
one additional assumption on Aq. 

We begin with the following observation: 

Lemma 3.2. Suppose that B C Q satisfies B C tp(B). Then B C Aq. 

Proof From B C <p(B ) it follows that ^ (B) C (/9 J+1 (H) for all j > 0. Hence 

B=f]P(B)c p| ^(Q) = A q . 


We now can prove our first result. 
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Proposition 2. Let Aq be the global attractor relative to the compact set Q , and 
suppose that the embedded attractor A & satisfies 4 cQ. Then 

(14) Ak C Aq. 

Proof. By construction of our dynamical system (see (5)-(7)), we have <p(A k ) = A k . 
Thus, Lemma 3.2 implies that A C Aq. □ 


Remark 3. Observe that we can in general not expect that Ak = Aq. In fact, by 
construction Aq may contain several invariant sets and related heteroclinic connec¬ 
tions. In this sense Ak will be embedded in Aq. 

Next observe that the Q^s define a nested sequence of compact sets, that is, 
Qe+i C Qi . Therefore, for each m, 

m 

(15) Qm = Pi Qu 

M 

and we may view 

oo 

(16) Qoo = f]Qe 

1=1 

as the limit of the Qi s. 

Now we will prove the convergence of the subdivision scheme for continuous ip. 
More precisely we will show that Qoq = Aq. This will be done in two steps. The 
first is 


Lemma 3.3. 


Q oo Aq. 


Proof. We will show that 

Q DC Cl p(Q oo)’ 

Then the result follows with Lemma 3.2. 

Let y G Qoo- Then for every i > 0 there is a unique Bi(y) G Bi with y G 
Bt(y). By the selection step of the subdivision scheme (see (13)), there is zi G 
Qi with p(zi) G Bi(y). Choosing a convergent subsequence of (z^), if neces¬ 
sary, we may assume that z = lim^oo By construction, z G Qoo, and since 
lirn^oo diam (Bi(y)) = 0 we conclude that lim^oo p{zf) = y. Finally p is continu¬ 
ous, and therefore y = p(z). Hence y G p(Qoo )• Cl 


For the inverse inclusion, we need to introduce an additional assumption, namely 
that p~ 1 (Aq) C Aq. This is automatically satisfied in the case where p is a 
homeomorphism. Moreover if Ak is attracting and Ak = Aq then Aq is backward 
invariant. These observations justify this assumption. 

Lemma 3.4. Suppose that p~ 1 (Aq) C Aq, then 

Aq d Q oo* 

This proof is identical to the proof of Lemma 3.2 in [5]. Thus, we will not restate 
it here. 

We summarize the convergence result in the following 

Proposition 3. Suppose that the relative global attractor Aq satisfies p~ 1 (Aq) C 
Aq. Then 

A-Q — Q oo- 
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3.4. Approximation of Attracting Sets. Observe that by construction of 9 ?, the 
set Aq defined in Section 3.2 contains the one-to-one image A k of the invariant set 
A of <f>. We now show that by using sufficiently high powers of <f> we can actually 
approximate a one-to-one image of A if A is attracting. 

Therefore, we now assume that A is an attracting set, that is, A attracts all 
bounded sets within a neighborhood U of A. Moreover we assume that the set Q 
is chosen in such a way that 

(17) A k C Q and E(Q) C U. 

Hence, for every x G Q, &(E(x)) will eventually approach the attracting set A for 
j -A 00 . However, this alone does not guarantee that A k is also an attracting set for 
the dynamical system ip. For instance, it may be the case that for a certain x G Q 
one has a “spurious fixed point” in the sense that 

x = (p(x) 

although &(E(x)) 7 ^ E(x) may be closer to A than E(x). 

In order to overcome this problem we now define for m > 1 the continuous maps 

(18) ip m = Ro$ rn oE 

and denote the corresponding relative global attractors by Aq. 

Remark 4. Observe that A is an invariant set for <F m for every m and therefore 
we can still use R as the restriction in our construction of the dynamical system 


Lemma 3.5. A k C Aq for all m > 1. 

Proof. Since <F(A) = A we have cp m (Ak) = A k for m > 1 . Moreover A k C Q (see 
(17)), and Lemma 3.2 implies that A k C Aq. □ 


Define 

n a q • 

m> 1 

Obviously A k C Aq. Moreover we have 

Proposition 4. A k = Aq . 

Proof. Suppose that x G Aq \A k - As A k is compact, this implies dist(x, A k ) = e > 
0. As A is compact, R is continuous and A k = ii(A), there is 8 > 0 such that 

dist(iz, A) < S dist(R(u), A k )) < 

Set V = E(Q) C U by assumption (see (17)). Since A is attracting and V is a 
compact set within U we can find an m > 1 such that 


h(4> m (H),A) <5, 


where h is the Hausdorff distance. By our choice of S it follows that 
h(R(<t> m (V)),A k ) = h(ip m (Q), Ak) < 


Thus, 


x fL Vm{Q) =* ® fL A% => X 0 


yielding a contradiction. 


□ 
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Remark 5. Roughly speaking Proposition 4 states that it will be possible to ap¬ 
proximate an attracting set for <f> if we perform the computations with appropriately 
high iterates of <f>. 

4. Numerical Realization for Delay Differential Equations 

As one important setting where finite dimensional dynamical phenomena occur 
in infinite dimensional Banach spaces, we consider delay differential equations with 
constant time delay r > 0. More precisely we consider equations of the form 

(19) y(t) = g(y(t),y(t - r)), 

where y(t) G M n and g : M n x M n -A M n is a smooth map. Following [14], we denote 
by C — C ([—t, 0], M n ) the (infinite dimensional) state space of the dynamical system 

(19) . Equipped with the maximum norm, C is a Banach space. 

Let y u (t ) be the trajectory generated by (19) with the initial condition u G C. 
Then the flow <f> s : C C of (19) is given by 

u i— y 4> s (r), where & s (u)(t) = y u (s — t) for t G [—r, 0]. 

Next we fix uj > 0 and consider the corresponding time-uj-map 4^ : C —>• C as our 
dynamical system. That is, we set 

(20) $ = and Y = C 
in our abstract dynamical system (1). 

In order to numerically realize the construction of the map (p = Ro<froE described 
in Section 3, we have to work on three tasks: the implementation of E, of R, and of 
respectively. For the latter we will rely on standard methods for forward time 
integration of DDEs [2]. The map R will be realized on the basis of Theorem 2.5 and 
Remark 1 by an appropriate choice of observables. For the numerical construction 
of the embedding E we will employ a bootstrapping method that re-uses results of 
previous computations. This way we will in particular guarantee that the identities 
in (6) are at least approximately satisfied. 

From now on we assume that upper bounds for both the box counting dimension 
d and the thickness exponent a are available. This allows us to fix k > 2(1 + a)d 
according to Theorem 2.5. 

4.1. Numerical Realization of R. For the definition of R we have to specify the 
time span uo and appropriate corresponding observables. In the case of a scalar 
equation (n * 1) we choose the observable / to be 

f(u) = u(-t). 

Thus, in this case the restriction R is simply given by 

R = D k [f, $](w) = (- u(-r ), $(w(-t)) } ..., $ fe_ 1 (w(-r))) T . 

The time span uj (see (20)) is defined to be a natural fraction of r, that is 

(21) co = L for K € N. 

K 

Remark 6. (a) Observe that a natural choice for K in (21) would be K = k — 1 

for fc > 1. That is, for each evaluation of R the observable would be applied 
to a function u : [—r, 0] M at k equally distributed time steps within the 

interval [—r, 0]. 
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(b) As described in Section 3.4 (see Remark 4) we will frequently replace by 
0 m (m > 1) in order to speed up the convergence towards the invariant 
sets A resp. A^. For an illustration of this procedure see Figure 2. 


O (U) 




Figure 2. Numerical realization of the restriction R for n = 1, K = 2 
and m — 6 K. 


For the numerical analysis of systems of delay differential equations (n > 1) we 
make use of Remark 1 as follows: For each component Uj of u we define a separate 
observable fj (j = 1 ,..., n) by 

(22) fj(u ) = Uj(vj) for a Vj E [—r, 0], 
and choose different time spans 

(23) Uj = -U for Kj e N 

accordingly. 

Finally, we note that also more general constructions for the restriction R : C 
R k can be employed. In fact, by virtue of Theorem 2.4, for any k that is sufficiently 
large for the delay embedding construction, an arbitrary linear map C -A will 
generically be one-to-one on A. Therefore almost every linear combination of tra¬ 
jectory points computed during forward integration can be used for the construction 
of the map R. 

4.2. Numerical Realization of E . In the application of the subdivision scheme 
for the computation of the relative global attractor A & described in Section 3.2 one 
has to perform the selection step 

B e = {b e Be : 3B e B t such that <p-\B ) 

(see (13)). Numerically this is realized as follows: At first cp is evaluated for a large 
number of test points z k E B for each box B Eft. Then a box B m is kept in the 
collection Bn if it is hit by (at least) one of the images p{z~). 

Remark 7. In practice the test points z k E Bk can be chosen according to several 
different strategies: In low dimensional problems one can choose them from a regular 
grid within each box B Alternatively one can select the test points from the 
boundaries of the boxes. In our computations we have sampled a fixed number of 
test points from each box at random with respect to a uniform distribution. 
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For the evaluation of(^ = Ro4>o£ T ata test point z we need to define the 
image E(z ), that is, we need to generate adequate initial conditions for the forward 
integration of the DDE (19). In the first step of the subdivision procedure, when no 
information on A is available, we proceed as follows. In the case of a scalar delay 
equation, that is n = 1, we construct a piecewise linear function u = E(z), where 

(24) u(U) = Zi, 

for ti = —t + i • cj, i = 0,. .., k — 1. Observe that by this choice of E and R the 
condition RoE(z) = z is satisfied for each test point z (see (6) and Remark 6 (a)). 

For n > 1 we proceed analogously and distribute the components of z e R k to 
the components Uj of u = E{z) G M n according to (22) and (23). Also in this case 
the condition R o E{z) = z still holds. 

In the following steps of the subdivision procedure we proceed as follows: Note 
that if B G S^, then, by the selection step, there must have been a Be Bi-\ such 
that R(^ UJ (E(z))) G B for at least one test point z G B. Therefore, we can use the 
information from the computation of <$> u (E(z)) to construct an appropriate E{z) 
for each test point z G B. 

More concretely, in every step of the subdivision procedure, for every set B G Bi 
we keep additional information about the points <& u (E(z)) that were mapped into 
B by R in the previous step. In the simplest case, we store /q > 1 additional 
equally distributed function values for each interval (— r + (i — 1)cj, —r + iuS) for 
i = l,...,fc — 1. When y>(B) is to be computed using test points from F>, we 
first use the points in B for which additional information is available and generate 
the corresponding initial value functions via spline interpolation. Note that the 
more information we store, the smaller the error || Q> u (E(z)) — E(z) || becomes for 
z = R(^ UJ (E(z))). That is, we enforce an approximation of the identity EoR{u) = u 
for all u G A (see (6)). 

If the additional information is available only for a few points in 5, we gener¬ 
ate new test points in B at random and construct corresponding trajectories by 
piecewise linear interpolation. 

5. Numerical Results 

In this section we present results of computations carried out for three different 
delay differential equations. In each case, u(t) is scalar, although for the DDE 
considered in Section 5.2 the problem is recast into a three-dimensional form in 
order to obtain a first-order equation. 

5.1. The Modified Wright Equation. As the first example, we consider a mod¬ 
ification of the Wright equation, 

(25) u(t) = — a • u(t — 1) • [1 — u 2 (t)\. 

In [14] it has been shown that the stationary solution uo(t) = 0 undergoes a super¬ 
critical Hopf bifurcation at a = tt/2. Thus, (25) possesses a stable periodic solution 
for a > 7r/2 - at least locally. In our computations we set a = 2, choose the em¬ 
bedding dimension k = 5, and approximate the relative global attractor Aq C M 5 
for Q = [—2,2] 5 , see Figure 3. Here the set Aq consists of a reconstruction of 
the two-dimensional unstable manifold of uq = 0 which accumulates on a stable 
periodic orbit at its boundary. 
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- 0.8 - 0.8 


Figure 3. Three-dimensional projection of an approximation of the rel¬ 
ative global attractor Aq within Q = [—2, 2] 5 for equation 
(25) after £ = 45 subdivision steps (embedding dimension 
k — 5 and iteration exponent m = 16, see Section 4.1) 


In Figure 4 we show box coverings of the reconstructed periodic solution itself. 
These have been obtained by removing a small open neighborhood U of the origin 
from Q = [— 2 , 2] 5 and computing Aq for Q = Q \ U . 

5.2. The Arneodo System with Delay. The second example is a modification 
of the Arneodo system [ 1 ] where a delay is introduced in the first order derivative 
of R, 

d 3 u . . d 2 u , . du , . , . o / \ 

^(f) + ^(f) + 2-(( - r) - au(t) +«(*) = 0. 

This equation has been introduced and analyzed in [21]. In our computations we 
use the equivalent reformulation as a first-order system 

u i = u 2 , 

= u 3 , 

ii 3 = —u 3 — 2u 2 (t — t) + aui — u\. 


(26) 
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(c) £ = 45 



(d) Simulation 


Figure 4. (a)-(c) Three-dimensional projections of successively finer 
coverings of the relative global attractor Aq which corre¬ 
sponds to a reconstruction of a periodic orbit of (25). (d) Pe¬ 
riodic orbit computed by direct simulation. 



(a) £ = 20 



- 0.8 - 0.8 
(b) £ = 30 


The undelayed system (i. e. (26) with r = 0) has been studied extensively. It 
possesses the equilibria Oi = (0,0, 0) and 0 2 = (cq 0, 0), the latter is asymptotically 
stable for a < 2. At a = 2 the equilibrium 0 2 undergoes a supercritical Hopf 
bifurcation (cf. [16]). For values of a which are slightly larger than two, points on 
the two-dimensional unstable manifold of 0 2 converge to the corresponding limit 
cycle on the branch of periodic solutions. That is, topologically speaking, we have 
the same situation as in Figure 3. 

For the delayed (i. e. r > 0) equation, the Hopf bifurcation occurs at decreasing 
values of a for increasing values of r. For fixed a = 2.5, the amplitude of the 
limit cycle grows with increasing values of r and loses its stability in a period¬ 
doubling bifurcation at r ~ 0.11 [21]. Our purpose is to investigate the structure 
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of the relative global attractor right after the occurrence of the period-doubling 
bifurcation. Concretely we set a = 2.5, r = 0.13, choose the embedding dimension 
k = 5, and approximate the relative global attractor Aq C M 5 for Q = [—1,5] x 
[-4, 2] x [-4,4] x [-4,4] x [—4,4]. This way we compute a reconstruction of the two- 
dimensional unstable manifold of the origin which accumulates on a period-doubled 
limit cycle. 

In this example we have made use of Remark 1 in our numerical realization. 
Concretely we have chosen uj = r /2 and the following three observables (see ( 22 ) 
and (23)) 

fi(u) = u 2 (-r), k\ = 3, A'l = 2, 
f 2 {u) = ui( 0 ), k 2 = 1 , 

/ 3 (m) = ^ 3 ( 0 ), k 3 = 1. 

Thus, the restriction R can be written as 

R(u) = (u 2 (-t),u 2 (-t/2),u 2 (0),u 1 (0),u 3 (0)) t 

Observe that R : C —> M 5 is linear and therefore also Theorem 2.4 could be used in 
order to justify this construction. The corresponding reconstructions of the relative 
global attractor are shown in Figure 5. 

Observe that after the period doubling bifurcation the relative global attractor 
contains a Moebius strip with the period-doubled periodic solution at its boundary. 
Thus, in the course of the period doubling bifurcation there has to occur a significant 
change of the geometry of the unstable manifold at its boundary so that it can 
accommodate the period-doubled solution. In fact, the corresponding mechanism 
has been analyzed analytically already in 1984 by Crawford and Omohundro [3]. It 
turns out that at the period doubling the unstable manifold starts to wrap itself 
“infinitesimally” around the unstable periodic solution. In a corresponding Poincare 
section this becomes a spiraling behavior with very sharp curvature, and we analyze 
this behavior at the reconstruction in Figure 6 (see also Figure 16 in [3]). However, 
we expect that one would have to choose a much higher resolution (i. e. higher 
number of subdivisions) in order to reveal this dynamical behavior more clearly. 

5.3. The Mackey-Glass Equation. Our final example is the well-known delay 
differential equation introduced by Mackey and Glass in 1977 [18], namely 

(27) ii(t) = B ^- 7 u(t), 

y J w M 1 + u(t-r) r > 1 

where we choose /? = 2 ,7 = 1, 77 = 9.65, and r = 2. This equation is a model of 
blood production, where u(t) represents the concentration of blood at time £, u(t) 
represents production at time t and u(t — r) is the concentration at an earlier time, 
when the request for more blood is made. Direct numerical simulations indicate 
that the dimension of the corresponding attracting set is approximately d = 2 . 
Thus, we choose the embedding dimension k = 7, and approximate the relative 
global attractor Aq for Q = [0,1.5] 7 C M 7 . In Figure 7 we show projections of 
the coverings obtained after £ = 28,42 and 63 subdivision steps as well as a direct 
simulation. 

Finally we conclude this section with an outlook and show a corresponding in¬ 
variant measure for the reconstructed Mackey-Glass attractor in Figure 8 . This 
measure has been computed with the software package GAIO [4] which is based 
on the techniques developed in [ 6 ]. However, a detailed investigation on how to 
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y 


/ x ^ (d) i = 45 with the period-doubled orbit 

' ' computed by direct simulation 

Figure 5. (a)-(d) Successively finer coverings of the relative global at¬ 
tractor for the Arneodo DDE (26) after subdivision steps 
(a = 2.5, r = 0.13, embedding dimension k — 5 and iteration 
exponent m — 15; see Section 4.1). 


approximate invariant measures in infinite dimensional problems efficiently using 
embedding theory will be done in future work. 

6. Conclusion 

In contrast to the situation for finite dimensional dynamical systems, for which 
there exists a wide range of advanced numerical tools, there are currently only few 
options besides direct simulation for the numerical computation of attractors of 
infinite dimensional dynamical systems generated e. g. by DDEs. In this paper we 
develop a general methodology for the computation of finite dimensional compact 
attractors of infinite dimensional dynamical systems, and illustrate its application 
to several DDEs. 

Combining the delay embedding technique with a set-oriented method for the 
computation of attractors, we obtain a flexible method for the analysis of infinite 
dimensional dynamical systems. More concretely, we use standard techniques for 
short-time simulation of the system in question to approximate the infinite dimen¬ 
sional dynamics, and employ the delay embedding technique on simulation results 
in order to obtain a representation of the dynamics by a continuous map on a 
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y x 

(a) Reconstructed unstable manifold. Red boxes 
(at y = 0) illustrate the location of the Poincare 
section. 



(b) Poincare section. Red circles mark the 
intersection with the period doubled periodic 
orbit. 



Figure 6. Detailed illustration of a Poincare section through the rela¬ 
tive global attractor for ( 26 ). 


moderately sized space. This map, in turn, can be analyzed using the subdivision 
scheme in order to compute a covering of an attractor. We show that in this way 
one obtains sets that are in one-to-one correspondence with the infinite dimensional 
attractor. 

The method proposed shares vital characteristics with its counterparts for finite 
dimensional systems. Most importantly, the numerical effort essentially depends on 
the dimension of the object to be computed, and not on the dimension of ambient 
space used for computations, that is in the case of this paper, the dimension of 
the space used for the delay embedding. In three examples, we have illustrated the 
suitability for the computation of two-dimensional attractors. (In the case of the 
Mackey-Glass equation, the attractor probably even has a box counting dimension 
larger than d = 2.6, see e.g. [10].) Furthermore, due to the set oriented nature of 
the underlying subdivision algorithm, the method does not depend on any special 
geometric properties (besides compactness) of the attractor. 
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